Anomalies in the Entanglement Properties of the Square Lattice Heisenberg Model 
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We compute the bipartite entanglement properties of the spin-half square-lattice Heisenberg model 
by a variety of numerical techniques that include valence bond quantum Monte Carlo (QMC), 
stochastic series expansion QMC, high temperature series expansions and zero temperature coupling 
constant expansions around the Ising limit. We find that the area law is always satisfied, but in 
addition to the entanglement entropy per unit boundary length, there are other terms that depend 
logarithmically on the subregion size, arising from broken symmetry in the bulk and from the 
existence of corners at the boundary. We find that the numerical results are anomalous in several 
ways. First, the bulk term arising from broken symmetry deviates from an exact calculation that can 
be done for a mean-field Neel state. Second, the corner logs do not agree with the known results for 
non-interacting Boson modes. And, third, even the finite temperature mutual information shows an 
anomalous behavior as T goes to zero, suggesting that T — > and L — >■ oo limits do not commute. 
These calculations show that entanglement entropy demonstrates a very rich behavior in d > 1, 
which deserves further attention. 



I. INTRODUCTION 

The study of bipartite entanglement properties in 
quantum statistical models is a promising way to un- 
derstand and classify their topological and universal 
propertiesii In one spatial dimension, the existence of 
finite entanglement entropy for gapped systems and a 
logarithmically divergent entanglement entropy for gap- 
less systems is well understood. The coefficient of the log 
divergence is universal in that it depends on the central 
charge and not on the microscopic details of the system. 
In two dimensions (2D), several gapless systems have 
been shown to display an area law, where the entropy 
per unit boundary length is necessarily non-universal, 
reflecting all the microscopic degrees of freedom at the 
boundary. Thus, looking for universal behavior necessi- 
tates looking for subleading terms. The area law coef- 
ficient will have a subleading power law dependence on 
the size of the system. However, more interestingly, there 
can be terms associated with corners on the boundary, 
with broken symmetry, or with non-trivial topology in 
the bulk which could be universal and allow one to clas- 
sify different topological and critical phases. 

The spin-half square-lattice Heisenberg model with 
nearest-neighbor interactions is one of the most stud- 
ied and best understood models in 2Di^ The model has 
long-range order, and spin rotational symmetry is spon- 
taneously broken in the ground state. This broken sym- 
metry is well described by a tower of rotor states;^ whose 
energy scales with system size L as l/L'^ in d dimen- 
sions. For large systems, these states are well separated 
from excitations around the ground state, which in the 
long- wavelength limit are spin- waves whose energies scale 
as 1/L. These spin-waves are believed to become non- 
interacting in the long- wavelength limit as long as one is 
away from a quantum critical point where the long-range 



order might go continuously to zero. 

In this paper, we have developed a number of different 
computational methods to calculate entanglement prop- 
erties, which are valid for arbitrary dimensional quan- 
tum statistical models. Stochastic series expansion (SSE) 
quantum Monte Carlo (QMC)^!^ and High Temperature 
Expansions (HTE) are methods that allow one to cal- 
culate thermodynamic properties of the model at finite 
temperatures, and can be used to obtain the Renyi mu- 
tual information associated with dividing the system into 
two regions A and B^ This mutual information should 
reduce to twice the entanglement entropy as T ^' 0. The 
Valence Bond (VB) QMC method is a particularly pow- 
erful tool for studying properties of the model directly at 
T — 0»i"— We have developed extensions of this method 
which employ loop updatesi^ necessary to accurately cal- 
culate the entanglement properties for finite lattices with 
different bipartite divisions. Finally, Ising expansions at 
T = provide yet another method to calculate proper- 
ties of the Heisenberg model in its ground state. By ap- 
proaching the ground state of the Heisenberg model from 
the ordered side one can calculate a specific ordered state 
and its entanglement properties. 

Our first result is an estimation of the leading non- 
universal area law coefficient, calculated as a = 0.097 ± 
0.001 in VB QMC (Section HIl) and a = 0.094 ± 0.001 in 
Ising series expansions (Section IIII[) . which are in good 
agreement. The small discrepancy shows that there are 
systematic errors not included in the numerical estimates 
of error bars. However, these are of order one percent. 

The universal contributions to the subleading scal- 
ing come from a number of sources. For example, one 
might expect the universal entanglement properties of 
the model to be related to the broken symmetry of the 
Neel state, and to the presence of free Bosons result- 
ing from non-interacting spin-waves. The entanglement 



properties of free Bosons (and free Fermions) have been 
computed in 2Dr^ Furthermore, the entanglement prop- 
erties of a mean- field Neel state, where all spins on sub- 
lattice 1 and sublattice 2 separately form a maximal spin 
state that are then combined into a singlet, can be calcu- 
lated exactly. This gives a mean-field bulk entanglement 
entropy associated with broken symmetry which scales as 
cln(€), with c = 2, where £ is the length of the boundary 
(Section IIV[) . In contrast, our VB QMC simulations for 
a bipartite division with no corners show a logarithmic 
term with c = 0.74 ± 0.02. Recently, the entanglement 
properties of this model were calculated using modified 
spin- wave theory;^ which gave an estimate of c = 0.92. 
Thus, our numerical simulations are much closer to spin- 
wave theory. This suggests that, in addition to broken 
symmetry, there are logarithmic contributions in the bulk 
that come from other sources. 

The corner contributions can be obtained in QMC by 
comparing the log terms in a system using a bound- 
ary with corners to a system without boundary corners, 
giving an estimate of —0.10 ± 0.02. The series esti- 
mate, —0.080 ± 0.008, gives reasonable agreement with 
the QMC results. In contrast, if one takes two free Bo- 
son modes contributing to the corners, then one gets 
- -0.0496 (SectionEl). 

Finally, we have studied the properties of the model at 
finite temperatures, and examined the approach to T — s- 
0. This can be done by SSE QMC and HTE. The HTE 
extrapolations for the Renyi mutual information agree 
well with QMC down to T « 1. Below this temperature, 
the QMC data shows a sudden decrease and a crossover 
to a lower saturation value as T — J- 0. The latter is 
consistent with the entanglement entropy calculated at 
T = 0. The HTE does not show any sharp decrease. The 
sharp decrease has a size dependence and could imply the 
limits of T ^^ and L — >■ oo do not commute. Such non- 
commutation is well known for other response functions 
of the Heisenberg modeU^i^^ and also suggests a sizeable 
non-mean-field contribution to the area law term. 



II. VALENCE BOND QUANTUM MONTE 
CARLO ALGORITHMS AND RESULTS 

In order to calculate the zero-temperature scaling of 
the Renyi entanglement entropy in the Heisenberg model 
on finite-size lattices, we employ the valence-bond (VB) 
quantum Monte Carlo (QMC) method developed by 
Sandvik.— This is a highly-efficient method to project 
out the model's ground state by repeated application 
of the Hamiltonian to a trial wavefunction, through a 
Monte Carlo sampling of bond operators. As an im- 
provement upon our previous entanglement measurement 
procedure^ we use a modified version of the more ef- 
ficient loop algorithm.— This modification allows for a 
change to occur in the Monte Carlo weight, by modify- 
ing the space-time topology of the simulation cell. As in 
Ref. [l5| , this modified weight is required so that one can 



measure the difference between entanglement entropies 
of two distinct system subdivisions, instead of the abso- 
lute entanglement entropy. Then, if the geometry of the 
regions are chosen properly the difference can converge 
faster than the bare entanglement entropy measurement. 



A. Valence Bond Quantum Monte Carlo 

First we briefiy discuss the basic VB QMC algorithm 
(for more detail see Refs. p-Q) and the more recently 
developed loop update (see Ref. [1^) which significantly 
improves the scaling of the algorithm. The foundation 
of the VB QMC technique is to project out the ground 
state of the system, done by applying a high power of 
the Hamiltonian H^^ to a trial state. We use the Heisen- 
berg Hamiltonian, rewritten in terms of bond operators 
[Hab = \ — ^a'^b) acting on pairs of sites (aand&), which 
are nearest-neighbor pairs in this paper. The Hamilto- 
nian to the power M can be written as a sum of possible 
arrangements of these bond operators in a list of size M . 
The Monte Carlo algorithm importance samples terms in 
this sum, using a weight which depends on the number 
of off-diagonal operators in the term.— 

In its original form, the VB QMC scheme can be 
used to project out one copy of the ground state wave- 
function (single-projector) or simultaneously project two 
copies (double-projector) which can be used to measure 
expectation values of observables in the simulationj^r— 
Our previous scheme for measuring Renyi entanglement 
entropyi^ employed a double-projector method for calcu- 
lating the expectation value of a Swap operator. In this 
paper, we develop a highly- efficient loop variation of this 
measurement algorithm, outlined below. 



1. The Loop Algorithm^ 

The loop update for VB QMC simulations was intro- 
duced in Ref. [1^ as a highly efficient way of carrying 
out the sampling procedure. In addition to working in 
a basis of valence bonds, this scheme also samples over 
spin states. This combined spin/bond basis is shown to 
eliminate the need for a rejection step, and thus samples 
operators and basis states with high efficiency. 

To begin, operators in this case are divided into two 
classes, 



Hab{l) = {\-SlSl) 

H,b{2) = -UStS^ + S-S^, 



(1) 
(2) 



called diagonal and off-diagonal operators respectively, 
where the sum i/ah(l) + Hab{2) is equal to the bond 
operators Hab from the standard VB QMC algorithms 
mentioned above.— ^— 

The loop algorithm is best visualized using a diagram 
of the simulation cell showing the placement of valence 
bonds, spins, and operators, as depicted in Figure[T] This 
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FIG. 1: A possible simulation cell diagram for a 6-site system, 
including loops, operators, the initial valence bond states, and 
the compatible initial spin states. Up (down) spins are shown 
in grey (black). (Off-) diagonal operators are shown in (blue) 
black, giving a total of M = 3 operators in this example. The 
dashed line denotes \Vl) (propagated from the trial state on 
left) and \Vr) (similarly from the right). 



trix. In contrast, it is typically a challenge to measure 
entanglement entropy with QMC, as the density matrix 
is not sampled in a straightforward way. Over the past 
few years there have been several proposals of entan- 
glement measures in Monte Carlo simulations^^"— Re- 
cently, methods to measure Renyi entanglement entropy 
using a swap (or more generally a permutation)^^ op- 
erator have been developed and implemented in several 
types of Monte Carlo methodsi^ii^ Next, we outline the 
basic methodology relevant for VB basis QMC,— based 
on the expectation value of the Swap operator, before de- 
tailing current advances for improving the measurement 
efficiency using a hybrid loop-ratio trick estimator. 



1. Renyi Entropies and the Swap Operator 



diagram represents two VB trial states, which are the left 
and right edges of the figure, each projected "inwards" 
(by M = 3 operators). The projected state occurs in 
the center of the diagram, denoted by the dashed line. 
Along with the trial VB states, initial spin states are se- 
lected at random, with the condition that the two spins 
in a single VB must be antiparallel. For each trial state, 
M operators are chosen such that they each act on a 
pair of antiparallel spins (in the initial step of the algo- 
rithm, these are all diagonal operators). There are then 
two types of updates in this algorithm: spin updates and 
operator updates. 

For the spin update, loops are first constructed by link- 
ing the operators and valence bonds (shown in Fig. [IJ. 
Then, for each loop that is built, a decision is made to 
flip all the spins in that loop, with probability 1/2. This 
update samples possible spin states for the given valence 
bond configuration. In the second type of update, the 
operators in the list are changed so that diagonal opera- 
tors are rc-samplcd at random, subject to the condition 
that they remain acting on antiparallel sites. This re- 
configures the propagated valence bond states and the 
topology of the simulation cell loops for future updates. 

Measurements can be computed as usual^ using the 
propagated valence bond states, \Vl) and |Vr), which 
can be extracted from the simulation cell diagram by 
following the loops crossing the dotted line in Fig. [TJ 

In the next section we discuss the measurement of 
Renyi entanglement entropies with either the double- 
projector or loop algorithms. In the following section 
we describe this measurement with the loop algorithm 
using a modified VB QMC simulation cell. 



B. Measuring Renyi Entropies with VB QMC 

Numerical techniques such as exact diagonalization or 
density matrix renormalization group (DMRG) simula- 
tions are able to directly measure entanglement entropy, 
since the calculation provides access to the density ma- 



Wc arc interested in the generalized Renyi entropies, 
which quantify entanglement between a system subdi- 
vided into two regions, A and B. They are defined as 
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■ln[Tr(pS)] 



(3) 



where p^ ~ Tr^ (p) is the reduced density matrix of the 
total system traced out over region B and the standard 
von Neumann entanglement entropy is recovered in the 
limit as a — ^ 1. 

Despite the inaccessibility of the full wavcfunction of 
the system in Monte Carlo techniques, it is possible to 
sample Tr(p^) for integer a > 1. This is accomplished 
by taking the expectation value of a Swap operator— for 
a = 2 or permutation operator for a > 2^ e.g. 



S2iA) 



-ln{{S'wapA)) 

-^in((n^)). 



(4) 
(5) 



To measure the a*"^ entropy, each projected state must 
be composed of a non-interacting copies of the system. 
The permutation operators 11^ act to cyclicly exchange 
the state in region A between the a different copies of the 
system, and are constructed such that (H^) = Tr(p^). 
In the case of spin states (which for each MC step in 
the simulation is a product state) we can simply swap 
the states within region A between copies of the system. 
For valence bond states the application of the permuta- 
tion operator also has a simple result; it acts to exchange 
the endpoints of valence bonds within region A between 
copies of the system, and as such can create bonds be- 
tween the non-interacting copiesJ^ 

The bare measurement of the Swap operator has been 
showi>i^ to have problems with convergence for large re- 
gion A, while in principle the measurement should be 
symmetric such that {Swap a) = {Swaps), since the two 
density matrices pA and pB have the same eigenvalues. 
This is because the exchange of a larger region gives a 
larger number of different states as a result of that swap. 



and thus a larger range and number of possible values in 
the expectation value. It simply takes more Monte Carlo 
steps to converge upon the same result. Another way to 
think of it is that, even though Tr(p^) = Tr(p^), if region 
A is much larger than region B then Tr(p^) contains ex- 
ponentially more terms than Tr(p^). Getting the value 
to converge by a stochastic sampling of these terms takes 
much longer. However, it is possible to instead sample a 
series of smaller regions and greatly improve the conver- 
gence time. This technique is described in the following 
section. 



2. The Ratio Tnck 

The convergence issue mentioned above can be ad- 
dressed by a reweighting of the Monte Carlo sampling 
scheme. Begin by considering the double-projector 
method; where one measures the expectation value of 
Swap by sampling terms from 

{SwapA) - ^ /T/IT/ (6) 

l^lrWlWr{Vl\Vr) 

where jV;), \Vr) are the states obtained by applying 
lists of bond operators to the trial states and wi , Wr 
are the weights accrued by applying those operators. 
Terms are sampled proportional to the total weight 
W = wiWriVilVr), by accepting a new configuration with 
probability W'^°'" /W°^'^. Thus one can simply measure 

{Vi\SwapA\V^[ 
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once per Monte Carlo step, and the average 

value will give us (SwapA)- 

The convergence difficulties can be combatted by using 
the ratio trickJ^ One modifies the sampling weight to 
include the expectation value of a Swap operator for a 
region A that is close in size to the region we intend to 
measure. One can then measure the ratio of these two 
operators, e.g. 

{SwapA) _ T.lr WrWr{Vl\SwapA'\Vr) lv!\Swapt!^v!) 



{Swap A') 



J2lr WiWr{Vi\SwapA'\Vr) 



(7) 
This improves the sampling since, if regions A and A' are 

similar in size, the measurement ,,, 'L'""^^ u^ \ will have 

' {Vi\Swapji,\Vr} 

fewer possible values than iv'lv) '' ' ^'^'^ those values 
will have a smaller variance. 

Note however that one is only measuring a ratio of ex- 
pectation values. That is, to obtain {Swap a), one must 
know the value of {SwapA')- If {SwapA') was also ob- 
tained by a ratio trick simulation, the expectation value 
for the smaller component of A' must be determined, and 
so on. Thus, the procedure that we use in this paper is 
to measure a range of sizes for region A, beginning with 
a measurement of the bare Swap for a small region size, 
and increase the size of regions A and A' over several 
simulations in sequence. That is, we measure Swap for a 
sequence of different region sizes, Ai, A2, ■ ■ ■ , A„, where 
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FIG. 2: One possible simulation cell configuration for the 
loop ratio algorithm on a 6-site system where a — 2 and region 
A' contains the first 2 sites of the system. Spins between the 
usually non-interacting copies are connected through loops 
via the Swap operator (red). The green and orange links are 
used to show the connections between the sheets. The loop on 
the left side of the swap operator on the top sheet is connected 
to the right side on the bottom sheet and so forth. 



the number of lattice sites in Ai^i is greater than the 
number of sites in region Ai . Then, the Renyi entropy of 
an arbitrary region A„ is calculated through 
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\{SwapA„_2) 

\n{{SwapA,)) (8) 



where each ratio is calculated via Eq. [71 and the last ex- 
pectation value for Ai via Eq. [51 Note that each term in 
the sum requires a different VB QMC simulation, since, 
although we can measure the entropy for any region A 
within one simulation, we can only use one size of A' per 
simulation, since it affects the sampling of the valence 
bond states as described below. The scaling cost of the 
Ratio trick is therefore n; however, the gain in sampling 
efficiency is demonstrated to more than compensate for 
this additional simulation cost. 



C. The Loop/Ratio Algorithm 

In order to calculate Renyi entanglement entropy with 
maximal efficiency in VB QMC simulations of the Heiscn- 
berg model, an algorithm should be employed that com- 
bines the loop update with the ratio trick. When modi- 
fying the loop algorithm to use the ratio trick, the same 
principles as in the double-projector algorithm (above) 
apply, however the sampling weight (from Eq. ([5])) is not 



explicit since one instead samples over spin states whose 
overlap is always unity. 

In order to made the necessary modification to the 
loop algorithm, the system should first be replicated 
so that two non-interacting copies are present, as usual 
for measurements of the Swap operatorj^ Then, links 
in the simulation cell are reconnected as if there were 
a Swap operator permanently applied to the projected 
state \Vr), shown in Fig.[2j This causes spins from differ- 
ent non-interacting copies of the system to be connected 
via loops, which means they can be flipped together, 
and thus the spin states are sampled according to the 
swapped system {Vi\SwapA'\Vr)- The measurement of 
{Vi\SwapA\Vr) / {Vi\Swapyi'\Vr) is then accomplished by 
measuring an operator which swaps the states of the sites 
in region A that were not already swapped in region A' , 
assuming A' C A. 

This method has the same limitation of the double- 
projector ratio trick, that only one value of A' can be 
used per simulation, so the region to be measured must 
be built up from a small region A according to Eq. ([8|). 
In our results below, we use two geometries for build- 
ing the region A: "strips" and "squares" (see Fig. U or 
Ref. 1201 ). Strips refer to geometries where region A has 
one dimension equal to the linear size of the toroidal sys- 
tem itself, and is therefore without corners. The Renyi 
entropy 52 (A) is built up through the ratio trick by sys- 
tematically adding sub-regions of size L x 1. Squares 
refer to geometries where the linear size of A is increased 
symmetrically, starting from size 1x1. Square regions A 
necessarily have four corners. 



D. Results 

We begin by testing some basic properties of the Renyi 
entropy as calculated through the VB basis QMC with 
the loop/ratio trick outlined above. First, we examine 
the convergence of the Renyi entropy as a function of 
operator list length per site, m — M/N, as illustrated in 
Fig. [21 Using several system sizes, both periodic and open 
boundaries, and strip and square geometries (described 
above and in Fig. |4]), we see very good convergence by 
m = 10. This value of m was used for all the following 
VB QMC simulations, as the measurements are able to 
converge with that number of operators, but additional 
operators would detract from the algorithm's efficiency. 

Figure m shows examples of 5*2 (A) for regions A of both 
strip and square geometries of increasing width for a 20 x 
20 toroidal system. The length of the boundary for all 
strip regions in this plot is £ = 40, whereas for the square 
regions £ = Ax. The area law scaling of 5*2 is apparent 
in that the strip and square geometries both approach a 
straight line with zero and non-zero slope respectively. 

To determine the scaling of entanglement entropy in 
two dimensions, we examine L x L systems with periodic 
boundary conditions, the results of which are shown in 
Fig.[SJ Region A was systematically built up according to 
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FIG. 3: Percent error in S2 versus the number of operators 
per site (m) for different L x L lattices. For L = 4, 6, 8 the 
lattices have open boundaries and region A is half the sys- 
tem using the strip geometry. The exact values were found 
using density matrix renormalization group (DMRG) simula- 
tions. The L — 20 lattice has periodic boundaries and A is 
a 2 X 2 square. The "exact" value is taken from the m = 10 
simulation. Each data set was fit to an exponential function. 



the square and strip geometries as defined in Fig. 31 Since 
each region is built to satisfy Eq. ^, one can perform 
fits for several sizes of subregion for each lattice size. In 
this case (as opposed to the plot in Fig. Uj) , for each set of 
data in Fig. [5] we use a region A with width proportional 
to the system size. This is done in attempt to overcome 
finite size effects, and the interaction of boundaries (that 
can be seen for the strip geometry in Fig. Sj). Fig. [51 
includes data for 5*2 {A) using regions with width x — L/2 
as well as smaller regions x < L/2. A smaller region A 
has the advantage that 5*2 (A) converges faster, but the 
drawback is it brings one into the regime of finite-size 
effects, apparent in Fig. [H with small x, where ^2 is lower 
than the area law value found at x = L/2. This deviation 
seems to depend on the fraction of the system contained 
within region A which is why, for the square case (where 
both the size of region A and the boundary length are 
scaled with system size), the different region widths do 
not change the entropy scaling very much. However, for 
the strip geometry with a region A width oi x = L/2 — n, 
the fraction of the system contained in region A changes 
with L, and approaches 1/2 as L increases, i.e. when 
L/2 ^ n. This effect is evident in the top panel of Fig. [SI 
where for n ^ Q the entropy is diminished for smaller 
system sizes, but approaches the n = values as system 
size increases. 

As is clear, the data gives excellent fits to the function 
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(9) 



where (. is the length of the boundary between regions A 
and B, provided that the very smallest lattice sizes are 
excluded. The values obtained for coefficients a, c, and d 
arc listed in Table [H 
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FIG. 4: (Left) The Renyi entropy for a 20 x 20 PBC sys- 
tem versus the width (x) of region A for both the square 
and strip geometries. The boundary length does not change 
with the region width for strip geometry (since its height y 
traverses the periodic lattice), thus the entropy approaches a 
constant value. For square geometry the boundary length is 
4a;. In both geometries the entropy scales with the length of 
the boundary, that scaling becoming better as the width of 
region A approaches half the system size. At right, the "strip" 
(top) and "square" (bottom) geometries on a 6 x 6 lattice. 
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FIG. 5: 82/1 vs £, where £ is the boundary length of region A, 
for regions with square and strip geometry embedded in Lx L 
systems with periodic boundary conditions. The data sets 
outlined in black correspond to region A with width x — L/2 
for both geometries. The other data sets use smaller region 
A of the same geometries. Fits to Eq. ((Qjl are included for all 
data, and the coefBcients found are listed in Table ID 



geometry 


L/2-X 


a 


c 


-d 


strip 





0.096(5) 


0.7(4) 


1.2(2) 




1 


0.095(6) 


0.7(9) 


1.3(6) 




2 


0.092(3) 


0.9(6) 


1.8(8) 




3 


0.08(8) 


1.2(0) 


2.(6) 




4 


0.08(2) 


1.5(4) 


3.(7) 


square 





0.097(6) 


0.64(2) 


1.0(6) 




1 


0.097(6) 


0.62(1) 


0.9(5) 




2 


0.097(7) 


0.61(7) 


0.9(1) 




3 


0.097(6) 


0.62(6) 


0.9(3) 




4 


0.097(5) 


0.6(3) 


0.9(4) 



TABLE I: The coefficients a, c, d found by fitting the data in 
Fig.[5]to Eq. (|9]). The fits were done for both strip and square 
geometries, beginning at a; = L/2 for both and decreasing the 
region sizes until x = L/2 — 4, where x is defined in Fig. U as 
the width of region A. 



III. SERIES EXPANSIONS 
A. Expansion Methods 

We have developed two different types of series ex- 
pansions to calculate the entanglement properties of the 
Heisenberg model. The first is High Temperature Ex- 
pansions (HTE) . This method was introduced in Ref. Q 
for the XXZ model, of which the Heisenberg model is a 
special case. In this method, the mutual information be- 
tween two regions A and B associated with their bound- 
aries or corners can be expanded in powers of inverse 
temperature /3. The calculation can be done for Renyi 
mutual information of index a by introducing a replicas 
of the system. The coefficient of /3" is a polynomial in 
a of order a — 1 so that the limit a — >■ 1 can be readily 
taken to calculate the von Neumann mutual information 
as well. 

Since there is no finite temperature phase transition 
in the Heisenberg model on the square lattice, the ex- 
pansions in /3 can, in principle, be extrapolated down 
to T = 0. It is well known that the correlation length 
of the system grows exponentially at low temperatures 
as exp(C/T). Hence, we expect corrections to lead- 
ing behavior to be exponentially small at low temper- 
ature, exp(— C/T). For this reason, a change of vari- 
ables w = tanh/3 is applied before Fade approximants 
are used. This extrapolation method has been used in 
the past for other properties of the Heisenberg model^i 
and also shows good convergence between different Fade 
approximants when applied to mutual information. 

In addition, we have developed a series expansion di- 
rectly for the entanglement entropy at T = using Ising 
anisotropy parameter A = Jxyl Jz of the XXZ model. 
When A = 0, we pick one of the two Ising states to ex- 
pand around. That unique state factorizes for any two 
regions A and B. So, any entanglement entropy vanishes 
in that limit. For any a > 2, the series expansion for 



the Rcnyi entropy can be calculated as a power series in 
A using a linked cluster expansioui ^^i^^ Unlike HTE, the 
coefRcients of the Renyi entropies for different a arc not 
related by a simple polynomial relation and for every a 
a different calculation is needed. Here, we will restrict 
ourselves to the calculation of Renyi entropy with a = 2. 
The Ising series expansions work in the thermodynamic 
limit, starting with a system that has a short correlation 
length and long-range order along a particular direction. 
Since a specific ordered state is picked out, their major 
limitation is that they cannot be used to study any bulk 
entanglement properties associated with broken symme- 
try. In any finite order of perturbation theory, contribu- 
tions can come only from the boundary between regions 
A and B. The series are non-singular as long as the sys- 
tem has a gap. The singular dependence on the size L 
of the system is replaced in the series expansion stud- 
ies by a dependence on the correlation length ^, which 
diverges as the gapless Heisenberg point is approached. 
By scaling, the dependence on L should translate into 
a similar dependence on ^. An advantage of series ex- 
pansions is that entanglement associated with different 
surface manifolds such as surfaces, lines and corners can 
be analytically separated and separate expansions can be 
obtained for the entropy associated with them. 



B. Results from Series Expansions 

Fade extrapolation of High Temperature Expansions 
(HTE) in the variable w — tanh/3 are compared with 
the data from stochastic series expansion (SSE) quantum 
Monte Carlo (QMC) simulations in Fig. [51 The agree- 
ment is excellent down to temperatures below J. At low 
temperatures the QMC data shows dramatic finite-size 
effects, with a maximum and a minimum, followed by a 
slow rise at lower temperatures. In contrast, the Fade ap- 
proximants show a steady monotonic rise and saturation 
at low temperatures. These results suggest that the limit 
of r — !■ and L — > cxd do not commute for these quan- 
tities. We note that such non-commuting limits are well 
known for other properties of the Heisenberg modeL ^'^i^^ 
However, these non-commuting limits have not been an- 
ticipated for the mutual information. In Fig. [71 we plot 
the crossover temperature, as defined by the temperature 
of the local maximum in mutual information at T > 0, as 
a function of L. We see that the crossover temperature 
scales inversely with log(X), or, equivalently, L scales ex- 
ponentially in the inverse temperature. Such a length 
scale agrees with the correlation length in the Heisen- 
berg model, which also scales exponentially at nonzero 
temperature, and so at first sight it is natural to ascribe 
this peak in the mutual information to the size of the 
correlation domains. One might imagine the following 
argument: when the size of the system is larger than 
the correlation domain size, there are many correlation 
domains along the boundary. Since the direction of the 
spin in a correlation domain can be viewed as roughly 



uncorrclatcd with that in other correlation domains, this 
contributes an additional term to the entropy of both re- 
gions. However, the correlation domains that cross the 
boundary produce some correlations between the two re- 
gions, leading to a positive contribution to the mutual 
information. While this argument captures the correct 
qualitative scaling, suggesting that the mutual informa- 
tion should drop at lower temperatures, it fails on quan- 
titative grounds. The contribution to the entropy of a 
correlation domain from the random ordering direction of 
that domain should scale something like the logarithm of 
the correlation volume (see, for example, the mean-field 
theory calculation of the next section, where at low tem- 
peratures the entire system comprises one correlation do- 
main), and hence should be proportional to /3. However, 
the density of these correlation domains (the number of 
domains per unit boundary length, which is inversely pro- 
portional to the correlation length) is exponentially small 
in f3, and so we should expect that this contribution to 
the entropy per unit length should become negligible as 
£ gets larger. In contrast, the numerical data shows the 
difference between the T > maximum and the T = 
limit increasing with increasing £, so that this difference 
remains as an unexplained phenomenon. 

Let the expansions around the Ising limit for the sec- 
ond Renyi entropy per unit length of the boundary be 
given by 

ri=2 

and, for a single corner, 

S2c = 2_^ &nA". 

The non-zero coefficients up to n = 14 are given in Ta- 
ble 2. Note that all odd order terms are zero, so the 
expansion is in the variable A^. The series for the line 
term iS'2/ are evaluated by first performing a change of 
variables S = X^ — 2A^ to remove a square root singu- 
larity at A = 1 and then calculating Fade approximants. 
From these we estimate 

52//^ = 0.094 ±0.001, 

for the Heisenberg model. Here, the error bars rep- 
resent spread in values obtained from different Fade 
approximantsi^ For the corner term, we expect a log- 
arithmic singularity of the form 



S2c = a::lnC 



-^ln(l-A2), 



where, we have used the fact that the correlation length 
diverges as (1 — A^)~^'^. Given the anticipated loga- 
rithmic singularity, we first differentiate the series with 
respect to the variable A^, and then use Fade approxi- 
mants biased at A^ = 1 to estimate the residue. From 
these, we estimate 

x = -0.020 ±0.002. 



n 


ffln 


bn 


2 


0.05555556 





4 


0.00314815 


0.00913580 


6 


0.00558342 


-0.00542753 


8 


0.00353554 


-0.00126847 


10 


0.00220329 


-0.00172570 


12 


0.00186784 


-0.00170418 


14 


0.00144690 


-0.00144091 



TABLE II: Ising series expansion coefficients for the line and 
corner terms 




QMC 8x8 
QMC 12x12 
QMC 16x16 
QMC 24x24 
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FIG. 6: The Finite temperature Renyi mutual information 
for L X L periodic systems with strip geometry (note £ — 2L) 
for region A, with SSE QMC resuhs for L = 8, 12, 16, 24, 32, 
and [4/6], [4/7], [5/5], [6/5] and [7/4] Fade approximants for 
the HTE in the variable w — tanh/3. 



Numerical values for S2i/i and x will be compared to VB 
QMC results and other calculations later in the discus- 
sion section. 



IV. MEAN-FIELD THEORY AND TOWER OF 
STATES MODES 

A. Tower of States and Thermodynamics 

One starting point for a theoretical treatment of the 
Hcisenberg model is a mean-field theory. This mean-field 
theory provides a simple framework for understanding 
some logarithmic bulk corrections associated with spon- 
taneous symmetry breaking, as well as providing some 
understanding of the low energy "tower of states" modes 
(first identified by Anderson^) and which are present even 
in the two dimensional model. 

In mean-field theory we consider the model Hamilto- 
nian 



H 



N 



7 . ^i ■ ^j ' 



(10) 
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FIG. 7: The inverse crossover temperature of Renyi mutual 
information for the system sizes shown in Fig. [6l where the 
crossover temperature is taken as the temperature at the high- 
est point of the T > mutual information peak. The a::-axis 
is logarithmic, showing that the crossover temperature scales 
as ~ log(I/). 



where the sum ranges over spin i in one sublattice and 
j in the other sublattice (we use 1 and 2 to denote sub- 
lattices), and there are total of N spins, with N/2 spins 
in each sublattice. The ground state of this Hamiltonian 
corresponds to taking all spins in sublattice 1 in a sym- 
metric state, with total spin N/A, and similarly taking 
all spins in sublattice 2 in a symmetric state, and then 
pair the two spins to form a singlet. 

We begin by working out thermodynamic properties of 
the model, and discuss their impact on QMC simulations 
using SSE and VB projector methods. We then turn to 
the question of the entanglement entropy. 

First, we work out the energy of the aforementioned 
singlet state as follows: let Si denote the total spin oper- 
ator in sublattice 1, and ^2 denote the same for sublattice 
2. Then, the Hamiltonian can be written as {J/N)Si ■82- 
This equals {J/2N)[{Si + 82? - Sf - 8^]. The singlet 
state has Si + 82 = 0. Looking at states with spin N/i in 
each sublattice, we have 8^ = 8^^ {N/i){N/4 + 1), so, 
the energy is -(J/iV)(iV/4)(A^/4 -I- 1), which is of order 
N. 

One can now look at excited states. The excited states 
where both sublattices have total spin N/4, but the whole 
system is not in a singlet, give the "tower states" : we find 
that in this case 5*1 4- 6*2 is not equal to zero. Indeed, the 
energy difference of this state, compared to the ground 
state, is simply given by 



E- E, 



grouU' 



d = {J/2N){8i + 82f 



(11) 



iei,je2 



That is, the energy is equal to {J/2N) times the total 
spin squared. In real two-dimensional systems, a similar 
low energy structure of states is observed, with the energy 
of these states proportional to total spin squared divided 
by N , although the coupling constant J describing the 
energy of these modes may be renormalized compared 



to the coupling constant appearing in the lattice Hamil- 
tonian. These low energy states can be observed in ex- 
act diagonalization, and in fact they are one of the best 
checks for the presence of symmetry breaking)^ Note 
that these states are much lower energy than the spin- 
wave states: in a 2D system with linear size L, they have 
energy of order 1/i^, while the lowest energy spin- wave 
has an energy of order 1/ L. 

There are also excited states where a given sublattice 
does not have total spin 7V/4. That is, not all spins in the 
same sublattice are in a symmetric state. One can check 
in this case that the energy of these states is increased 
above the ground state energy by an amount that is of 
order 1 (or more, if the spin is reduced a lot compared to 
iV/4) . These states can be viewed as a mean- field theory 
analogue of the spin-wave states; that is, the mean-field 
theory raises the energy of the spin- waves from order 1/ L 
to order 1. 

So, at a temperature of order 1, it becomes reasonable 
to ignore those states in mean-field theory (more precisely 
we need a temp of order 1/ In(A^)). In a two-dimensional 
system, the temperature needs to become of order 1/L 
to ignore the spin-wave modes. 

We now consider the effect of the tower of states on 
the bulk entropy, at temperature sufficiently low that the 
spin-waves can be ignored. We emphasize that this es- 
timate is not a calculation of the entanglement entropy, 
but rather a calculation of a bulk thermodynamic en- 
tropy. Considering the tower states, the total entropy is 
not hard to work out: the allowed states are 1 state (the 
ground state) which is a singlet, 3 states with spin 1 and 
energy increased by JS[S + \)/2N = 2J/2N, 5 states 
with spin 2 and energy 6J/2N, etc... (in this case, the 
Clebsch-Gordon coefficients work out simply, so that the 
number of states with spin S is exactly 25" -I- I). As a 
rough approximation, at a temperature T, we expect to 
excite states with spin such that JS'^/N is of order T. 
So, 5' is of order y^TN/J or less. The number of such 
states is of order S'^, and hence of order TN/J. Thus, at 
a temperature T much larger than 1/N, but sufficiently 
small that the spin-waves can be ignored, the entropy is 
equal to 



log(TiV/ J) + const. 



(12) 



B. Effect of Tower of States on QMC Simulations 

In QMC simulations, if we want to access the ground 
state, it is important to access a temperature sufficiently 
small that these tower of states modes can be ignored. 
Since the energy of the tower of states modes is so small 
(of order 1/N) this can require a prohibitively small tem- 
perature. However, we will see that this is not as big a 
problem as it might seem. 

First, consider the VB projector method. The pro- 
jector method starts with a trial wavefunction and then 
applies a large power of the Hamiltonian to this wave- 



function. The Hamiltonian is appropriately scaled such 
that its ground state will have the largest eigenvalue, so 
this high power of the Hamiltonian, acting on the trial 
wavefunction, produces a state close to the ground state. 
This is essentially the power method of finding largest or 
smallest eigenvalues of a matrix. Given that the Hamil- 
tonian has tower of states modes with energy very close 
to the ground state, it would seem that we would have 
to apply a very high power of the Hamiltonian, a power 
which is of order N, in order to produce a final state 
close to the ground state. However, the trial wavefunc- 
tion in the VB projector method has total spin 0, and 
the Hamiltonian conserves spin. Thus, the wavefunction 
produced by acting on the trial wavefunction with a high 
power of the Hamiltonian also has total spin 0. How- 
ever, all the tower of states modes have non-zero spin, 
and hence the wavefunction wc produce has no overlap 
with the tower of states modes. This is the reason why it 
suffices to simply go to high enough order to project out 
the spin-wave modes, as once the spin-waves excitations 
are projected out, the tower of states are also projected 
out by symmetry (note that we do have excited states 
with total spin which excite both a spin-wave and a 
tower mode, but such states have energy of order 1/L or 
higher, not 1/iV). 

Consider finally the SSE method. As discussed above, 
the tower of states modes do have a noticeable effect on 
the calculation of the bulk entropy. However, as we will 
see below, at least in mean-field theory they have only 
a small effect on the calculation of the entropy of the 
reduced density matrix, suggesting that the SSE calcu- 
lations of the reduced density matrix entropy converge 
well even without accessing temperature of order 1/A'^. 



C. Entanglement Entropy 

We now consider the entanglement entropy in the 
mean-field model. The entanglement between the 1 and 
2 sublattices is the simplest to calculate. Each sublattice 
has total spin N/A, and hence has A^/2 -I- 1 states. The 
ground state, which we call ip'^, is maximally entangled, 
and hence has entanglement entropy ln(A^/2 + 1), with 
all different Renyi entropies equal. 

The more interesting entanglement entropy to calcu- 
late is to imagine dividing the system into two halves, A 
and B, with each half having N/A spins in sublattice 1 
and N/A spins in sublattice 2. This is, in our opinion, 
the simplest mean-field model which is relevant to the nu- 
merical calculations on the two dimensional Heisenberg 
model, as in that case each region A and B contains spins 
from both sublattices. We will see that this calculation 
gives a logarithmic dependence on N also, although the 
result is more complicated. 

The ground state ip'^ can be obtained by taking any 
state which has spin N/A in each sublattice and which 
has non-zero overlap with the ground state and project- 
ing it into the total spin zero sector. We choose to use a 
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Nccl state as our state before projection, where we define 
this state ip^°°^ as a state in which ah spins in sublattice 1 
are pointing up and aU spins in sublattice 2 are pointing 
down. We do this projection by averaging over different 
rotations of the Neel state by all possible rotation matri- 
ces in the group SU{2). This may be formally written 



Ju£SU{2} 



KW iNcol 



(13) 



where the integral is over all unitaries U in SU{2) with 
the Haar measure and f/*^^ denotes acting with the same 
rotation U on each spin iV, and Z is a normalization 
factor. Writing -ip^ in this form makes it explicit that it 
has spin since all states with non-vanishing spin are 
projected out by the rotations. 

Since the Neel state we choose is oriented along the z- 
axis, in fact we can reduce the number of parameters in 
the integral from three (the number of angles needed to 
specify an SU{2) rotation) to two (the number of angles 
needed to specify the direction of a single spin), and so 
we write 



^0^^-1/2 / d6'd^i?(6', (/))?A 



Noel 



(14) 



where i?(0, (/)) is eqaal to exp{iJ2j{Sai+(f>a0) (i.e., it ro- 
tates each spin j about x and y axes and the sum ranges 
over all spins j) and the measure is chosen so that we ob- 
tain a uniform superposition of all possible orientations. 
The particular parameterization of rotations is not very 
important since we will be interested in the case that 
angles are small, discussed below. 

We compute the second Renyi entropy, 5*2. The calcu- 
lation of other Renyi entropies are similar; in this case, in 
contrast to the previous entropy calculation above, differ- 
ent Renyi entropies differ, so that 5*2 is not equal to the 
von Neumann entropy. To calculate the Renyi entropy, 
we must calculate 



where Swap a is the swap operator used in Ref. [if 
This expectation value is equal to 

Z-2 /d6'id(/)id6i2d<?!)2d6i3d03d6'4d(/)4 



(15) 



(16) 



X (V-N^^i ® i^N*°i|i?(0i, </)i)t ® R{e2, q^2y\SwapA\ 

i?(03>3) ® i?(^4, 04)1^''''"' ® ^'"''='). 

We first estimate Z as follows. We have 



(17) 

The integral / d6'2d</-2(V'^'^°1i?(6'i, 0i)^i?(6l2, </'2)|V''~''^''') is 
independent of 6*1,01. This can be seen as follows: 
the state J d62d(f>2R{02,(l>2)\'4'^^^'') has spin 0, so it has 



vanishing inner product with any state that has non- 
vanishing spin and hence it has the same inner product 
with the unprojected state (■!/;'^''°'| as it does with the 
projected state /d0id(^i('0^''°'|i?(0i,0i). So we can fix 
6»i = 01 = giving 



z= d(t)d0{ip^'"''\R{e,(t>)\ij^'"'') 



(18) 



The expectation value in the above integral is just the 
A^-th power of of (t|i?(0, 0)|t). This is approximated 
by (for smah 6*, 0) cxp{-N{e^ + (j)^)/2). Because of the 
factor of N in the exponent, the restriction to small 9, </> 
is justified as the expectation value is negligible for large 
6*, (j). Then, the integral over 0, is Gaussian and the 
result is that 



ZociV' 



(19) 



Now we estimate the integral in Eq. ((T5)) . The evalua- 
tion is similar, but in this case, we have an integral over 
four pairs of angles. Again, we can fix one pair of angles, 
such as 9i,(j)i, to equal 0. The overlap is a product of 
two overlaps, one for spins not in region A and one for 
spins in region A. The overlap of spins not in region A 
is small unless 9i is close to 6*3 and 0i is close to 03 and 
also 6*2 is close to ^4 and 02 is close to 04. Similarly, 
for the overlap of spins in region A to be non-negligible 
we need 9i close to 9^ (and 0i close to 04) and ^2 close 
to 9 3 (and 02 close to 03). Thus, the overlap of spins 
in region A forces certain pairs of angles to be close and 
the overlap of spins not in region A forces other pairs of 
angles to be close. Combining these conditions we see 
that all the angles need to be very similar. Then, we get 
(approximately) a Gaussian integral over three pairs of 
relative angles. The result is proportional to 1/N^, then, 
up to constant factors. 

Thus, the expectation is l/N^/{l/N^) = 1/A^, giving 
again an entropy which scales as In(A^) -t- const, for large 
N. Finally, we can consider this entropy at a non-zero 
temperature. The effect of a non-zero temperature (high 
enough to excite the tower modes but low enough to avoid 
exciting spin-wave states so that the total spin in each 
sublattice will still equal A^/4) is to give a global density 
matrix 



p = z- 



d0ld0id02d02F(01,01,02,02) (20) 



Xi?(01,0l)|^N-')(V^N-l|i?(02,02)n 

for some function F, which depends only on the rela- 
tive angles between the two rotations. That is, at higher 
temperatures the rotation angles in the bra and ket vec- 
tors become coupled, while at zero temperature F is a, 
constant. One can check that this still leaves us with 
an entropy of region A which is equal ln(7V) -|- const., so 
that, at least in mean-field theory, the In(iV) term in the 
entropy of region A does not depend upon whether or 
not the tower modes are excited. 
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V. SPIN- WAVE THEORY 

In understanding the area law for the Heisenberg 
model, an immediate question arises: what is the en- 
tanglement entropy in spin-wave theory? After all, spin- 
waves are gaplcss, and gapless modes might make one 
concerned whether or not an area law holds. In fact, the 
entanglement entropy of gaplcss modes depends strongly 
upon dimension. Free bosons with a linear disper- 
sion relation have a logarithmically divergent entangle- 
ment entropy in one dimension, following conformal field 
theoryji2^ but in two or more dimensions they obey an 
area lawi^ 

In Ref. [23, it was shown numerically that the Heisen- 
berg model itself obeyed an area law using density matrix 
renormalization group and a spin-wave calculation also 
led to an area law. In Ref. Il2l . a spin- wave calculation 
was carried out for a finite size system and was shown to 
roughly match the qualitative behavior from a quantum 
Monte Carlo simulationj^^ 

While gapless bosons with linear dispersion do pro- 
duce an area law in two dimensions, they also produce 
nontrivial exponents associated with corners<^ The en- 
tanglement entropy of a region A has a term equal to the 
log of the length scale of A, multiplied by the sum over 
corners of a scaling function of angle of each corner. For 
the entropy 6*2, this term is equal to 



-0.0062 ln(£) 



(21) 



for each 90-degree corner for a real scalar field with lin- 
ear dispersioui^ Note that the sign of this correction is 
negative. 

For the system wc are concerned with, we must multi- 
ply this result by two. There are two ways to understand 
this counting of modes, to see why the result must be 
multiplied by two. On the one hand, we can consider an 
0(3) nonlinear sigma model. In 2 4- 1 dimensions, this 
model has a symmetry broken phase, and the Heisenberg 
model ground state corresponds to this phase. In the 
symmetry broken phase, there are two Goldstone modes, 
corresponding to two different transverse directions in 
which the order parameter can move. In the Hamilto- 
nian spin- wave language, this factor of two again arises, 
but for a reason which initially might seem to be differ- 
ent. Suppose wc use a spin- wave representation in which 
the operator bj always creates an excitation on site i (we 
arc choosing to follow the notation of Ref. [ij, though of 
course the Hamiltonian spin-wave calculation is a text- 
book calculation); that is, if wc do a spin wave expansion 
about a state with spins up on the 1 sublattice and down 
on the 2 sublattice. then this operator bj corresponds to 
a lowering operator on the 1 sublattice and a raising op- 
erator on the 2 sublattice. Then, we find gapless modes 
at momenta near (0, 0) and (tt, tt) on a two dimensional 
square lattice. Thus, we again see a factor of two, arising 
from the existence of two different gapless points. 

We can clarify the relation between the factor of two 
in these two different approaches. These gaplcss modes 



in the Hamiltonian model correspond to the following 
states, respectively. Acting with the operator J2i ^I on 
the ground state, which creates a zero energy excitation 
with momentum (0,0), produces a state which is a super- 
position of all possible ways of flipping one spin. Acting 
with the operator X]i(~l)*^I' which creates a zero en- 
ergy excitation with momentum (tTjTt), produces a state 
which again is a superposition of all possible ways of flip- 
ping one spin, but with a plus sign for flipping a spin in 
the 1 sublattice and a minus sign for flipping a spin in the 
2 sublattice. The flrst of these states corresponds to act- 
ing with the operator ^- af on the ground state and the 
second to acting with the operator ^- af on the ground 
state. Thus, they correspond to two different ways of ro- 
tating the symmetry broken ground state, either in the 
YZ plane or the XZ plane, matching the two different 
Goldstone modes above. 

Therefore, for a square region, which has 4 such cor- 
ners, wc expect a correction of 



-0.0496 In(^). 



(22) 



As we saw above, numerical results disagree with this, 
suggesting some nontrivial effects that are not accounted 
for by this framework. 



VI. SUMMARY AND DISCUSSION 

In this paper, we have presented several computational 
methods to calculate the entanglement properties of lat- 
tice statistical models, in dimensionality greater than 
one, in a systematic manner. The stochastic series expan- 
sion QMC and high temperature expansions are finite- 
temperature methods. Since the former works with fi- 
nite systems, it is possible to calculate the ground state 
properties by going to sufficiently low temperatures. The 
latter is a series expansion, defined in the thermodynamic 
limit. It can, in principle, be extrapolated to T — >■ limit 
if there is no finite temperature phase transition. An in- 
teresting finding of our paper is that the limits T — > 
and L ^f OQ need not commute in these calculations. 

We have also developed computational methods that 
work directly at T = 0. The valence bond QMC is an 
exact ground state projection method for a finite system. 
The Ising series expansions represent an expansion in ex- 
change anisotropy around a given classical state. Using 
these methods we have calculated the area law associ- 
ated with the boundary, as well as sublcading logarithmic 
terms associated with the corners and the bulk. When- 
ever quantities are calculated by two different methods, 
there is good quantitative agreement. 

We have also discussed a mean-field calculation of 
the entanglement properties as well as the results ex- 
pected from non-interacting Bosons. Here, we find a 
surprise that the numerical results disagree with sim- 
ple expectations. The mean-field state, where all spins 
on one sublattice are equally entangled with spins on 
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the other sublattice, predicts a bulk log term of 2\ni, 
whereas our valence bond QMC results give cln£, with 
c = 0.74 ± 0.02. The latter is closer to a recent spin- 
wave calculation, where an estimate of c = 0.92 was 
obtainedf^ The discrepancy is even more puzzling for the 
corner terms. If long- wavelength spin- waves act as non- 
interacting Bosons, they should contribute a log term 
with a coefficient of —0.0496. That number should be 
compared with —0.080 ± 0.008 obtained in Ising series 
expansions and « —0.1 in the QMC studies. These dis- 
agreements suggest the need for further theoretical study. 
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